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Abstract 

We introduce and investigate a series of models for an infection of a 
diplodiploid host species by the bacterial endosymbiont Wolbachia. The 
continuous models are characterized by partial vertical transmission, cy- 
toplasmic incompatibility and fitness costs associated with the infection. 
A particular aspect of interest is competitions between mutually incom- 
patible strains. We further introduce an age-structured model that takes 
into account different fertility and mortality rates at different stages of 
the life cycle of the individuals. With only a few parameters, the ordinary 
differential equation models exhibit already interesting dynamics and can 
be used to predict criteria under which a strain of bacteria is able to in- 
vade a population. Interestingly, but not surprisingly, the age-structured 
model shows significant differences concerning the existence and stability 
of equilibrium solutions compared to the unstructured model. 

Keywords Wolbachia, cytoplasmic incompatibility, age-structured 
population dynamics, stability analysis 

1 Introduction 

Wolbachia is a maternally transmitted bacterium that lives in symbiosis with 
many arthropod species and some filarial nematodes [HI [2S] . It inhabits testes 
and ovaries of its hosts and has the ability to interfere with their reproduc- 
tive mechanisms, resulting in a variety of phenotypes. Well known effects are 
cytoplasmic incompatibility, induction of parthenogenesis, and feminization of 
genetic males, depending on the host species and the Wolbachia type. Besides 
the intrinsic interest in these mechanisms, Wolbachia are investigated as tools 
to drive desirable genes into a target population [T7] , as reinforcers of speciation 
[T2j [20l \21\ , and as potential means of biological control [14] . It was recently 
shown by McMeniman et al. [M] that infection with Wolbachia shortens the 
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life span of the mosquito Aedes aegypti, a vector for the Dengue fever virus. 
Since only older mosquitoes are carriers, this is a promising strategy to reduce 
the transmission of pathogens, without the ethically untenable eradication of a 
vector species. 

Beginning already a half a century ago [T] , various mathematical models for 
the spread of a Wolbachia infection have been proposed and studied in the liter- 
ature, see e.g. [lEl[Il|I71[iai2nilMI2S]and references therein. Largely, these 
models fall into two classes, depending on whether time proceeds in discrete 
steps or continuously. Examples for continuous models are the papers [12] and 
|18j that employ ordinary, respectively partial differential equations (with a spa- 
tial structure in the latter case). In their paper [T2], the authors proposed and 
studied a simple continuous model for the infection of an arthropod population 
with cytoplasmic incompatibility (CI) causing Wolbachia. Cytoplasmic incom- 
patibility in diplodiploid (i.e. with diploid males and females) species manifests 
itself in completely or partially unviable crosses of infected males with unin- 
fected females. For a discussion of the more complex outcome of cytoplasmic 
incompatibility in haplodiploid species see [53]. 

In this paper, we introduce a series of models for different aspects of inter- 
est. We start in section [2] with an ordinary differential equation model for a 
single Wolbachia strain that infects a population without separate sexes. In 
section [3] we present a model for infections with multiple strains. The present 
theoretical literature offers a complex picture of infection with multiple strains. 
While some authors exclude the coexistence of multiple strains of Wolbachia in 
infected individuals [Q] , others model doubly infected individuals as a class 
of their own plfSSl. Moreover, different assumptions are made about the mutual 
compatibility of individuals carrying different strains. We construct a general 
model that encompasses these different possibilities by suitable choices of pa- 
rameter values and/or initial conditions. Finally, motivated by the study |14j . 
in section [4] we refine our model from section [2] by considering age-structured 
populations. We refer the interested reader for basic concepts and results in 
structured population dynamics to [S] [T21 |M]- Modeling structured popula- 
tions usually involves partial differential equations which are more difficult to 
analyze. Analytical progress is still possible, and as we will see in [4] the age- 
structured model may exhibit richer dynamics. At this point it will be possible 
to study age-dependent killing of Wolbachia infected individuals. Our models 
contain parameters of only three types, namely transmission efficacies, levels of 
cytoplasmic incompatibility and fitness costs for the infected individuals. The 
analysis of the models aims to give conditions for the stability of specific equi- 
librium solutions that correspond to successful invasions. The paper ends with 
a discussion in section [5] and an outlook on future work. 
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2 Single-sex model for a singular Wolbachia 
strain 



Assume that the ratio of infected males to infected females is the same as 
the ratio of uninfected males to uninfected females and hence the population 
can be formally considered hermaphroditic. Let / and U denote the number of 
infected, respectively uninfected individuals in the population. Vertical trans- 
mission is partial, let r € [0, 1] be the fraction of infected offspring from infected 
parents (another common notation is fj, = 1 — r for the fraction of uninfected ova 
produced by an infected female, see e.g. [El |23]). Furthermore, we follow 
Keeling et al. |12| and assume that the birth rate for both infected and unin- 
fected individuals is equal (no reduction in fecundity in infected individuals). 
Let this rate be denoted by 6 > 0. Death of the individuals is modelled by a 
logistic loss term with rate d > that accounts for competition among the total 
population. However, infected individuals can suffer an additional loss of fit- 
ness given by _D > 0. Cytoplasmic incompatibility arises when an infected male 
fertilizes an egg from an uninfected female. Then, with a probability q £ [0, 1], 
the offspring is nonviable. As we do not consider separate sexes in this simple 
model, we just reduce the amount of offspring from uninfected individuals based 
on the probability of an encounter with an infected individual. Uninfected indi- 
viduals still arise due to incomplete transmission of the bacterium from infected 
parents. Our equations read 



- d{I + U)j U. 
-- , we obtain the reduced 

(1) 

{i + u)^ u. (2) 

Observe that t]~^ =: £, £ (0, 1] and that 1 — ^ can be interpreted as the fitness 
cost associated with Wolbachia infection. The point (0, 0) can be added to the 
domain of the state space, with the understanding that it is an equilibrium 
solution. It is obvious that the subspace of uninfected populations {0} x K 
is forward invariant (that is, if initially there are no infected individuals, then 
there will be none at later times) and if transmission is complete, t — 1, then 
so is the subspace of completely infected populations M x {0}. 
Equation ([T])-([2]) always admits the disease free equilibrium 

(io,uo) = (0,1). (3) 



- = (rfe-(d + /?)(/ + [/))/, 



Upon rescaling the time by t — > 5t and setting rj 
system for the quantities i = dl/b, u = dU/b 

^ = {t - T]{i + u))i, 

' i 



— = 1 - r)^ + 1 - g, 
dt ^ ' V ^ + ^t 
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Setting the left-hand side of equation ([T]) to zero and solving for an equilibrium 
point u yields 

u = — i. (4) 

Inserting this into the equilibrium condition for ([2| gives a quadratic equation 
for I, 

^i^ + iTi(-l)-q)l + T^il-^T) = 0. (5) 

Provided that 

(r(C-l)-g)2-4g(l-eT) >0, (6) 
equation ([s]) has the solutions 

< (q - r(e - 1) - V(^(^-l)-'Z)'-4g(l-eT)) 

H = Hiq,T,0 = ^ , 

(7) 

re (q ~ r(e - 1) + ^ (T{i - \) - q^ - Aq{\ - ir)) 
«2 = i2('7, r, = , 

these are always non-negative. The corresponding equilibrium values for the 
uninfected individuals are given by Q. For it2 > it is necessary that 

g + r(e - 1) > 0. (8) 

This condition is also sufficient, since then one can derive from r < 1 the 
inequality 

q + r(e - 1) > ^/{T{i - 1) - q)2 - 4g(l - ^r) 

and hence > 0. Condition ^ separates two regions in the parameter space, 
depending on whether other equilibrium solutions than the disease free equilib- 
rium are possible (see figure [T]) . We calculate the Jacobian of the right hand 
side F of Q-llI), 

CT — ri{2i + u) —rji \ 



At the disease free equilibrium (|3]), we have 

DFiO,l) = 



(9) 



T — T] 

— T — q —1 

this matrix has the eigenvalues —1 and t — rj < 0. The latter eigenvalue is 
only if r = 1 (complete transmission) and D = (no penalty for infection), in 
all other cases it is strictly negative, and the disease free equilibrium is locally 
asymptotically stable. 

Explicit expressions (with respect to the parameters q, r and f ) can be 
obtained for the eigenvalues of the Jacobian Q at the equilibrium solutions 
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Figure 1: The yellow surface separates the t, q)-paranieter space of model 
([T])-([2]) into a region where only the disease free equilibrium ^ exists (A) and 
where coexistence equilibrium solutions and (u2,«2) given by ([t]) are 

possible (B). However, only in the region (C) above the blue surface given by 
^ is U2 — tS^ — i2 > (this belongs to the observed stable equilibrium {u2,i2))- 



using e.g. mathematica (the mathematica notebook is available from the 
authors upon request). Unfortunately these expressions are very complicated 
and not easily analyzed. We will present instead some representative numerical 
examples to demonstrate the possible behaviours. 

Example 2.1. Assume that the infection is completely inherited, r = 1, cy- 
toplasmic incompatibility is complete (7 = 1, and that the cost of the infec- 
tion is low, ^ — 0.9. Then the three equilibrium solutions (io,ito) = (0,1), 
(ji,Mi) — (0.09,0.81) and {12, U2) — (0.9,0) are present, of which (jqiWo) and 
(12, M2) are locally asymptotically stable. The vector field is shown in figure [2] 
(left panel). The epidemic equilibrium (12, M2) has a much bigger basin of at- 
traction than the disease free equilibrium (iq, uq), in other words, the threshold 
for an infection to establish itself is low. 

Example 2.2. At high levels of cytoplasmic incompatibility, q — 1, and no 
penalty for the infection, ^ = 1, and a partial transmission r — 0.76, besides 
the equilibrium (io,Mo) there is the locally asymptotically stable coexistence 
equilibrium (12,^2) = (0.456,0.304), see figure |2] (right panel). The basin of 
attraction of the disease free equilibrium (ig, uq) is much bigger than in example 
|2.1[ indicating that the threshold for establishment of the infection is higher. 

Example 2.3. If the cytoplasmic incompatibility is very weak, q — 0.1 and 
the fitness cost of the infection is higher, ^ = 0.5 (and again r = 1), then 
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Figure 2: (Left) The vector field ([T])-([2]) for the parameter triple {£,,T,q) = 
(0.9, 1, 1) together with the three equilibrium solutions. Solid disks indicate 
locally asymptotically stable equilibrium solutions, while the disk indicates an 
unstable equilibrium. Shown are also regions of growth of u (light blue) and 
growth of i (light red). The blue lines are the stable manifolds of the saddle point 
(zi,Wi) and the separatrices of the equilibrium solutions («o,uo) and («2jW2)- 
(Right) The vector field for the parameter triple (^, r, g) = (1,0.76,1), 

which admits bistability and true coexistence. The stable manifold of the saddle 
point {ii, ui) is shown in blue. 



the equilibrium (io,uo) — (0,1) is globally asymptotically stable in the set 
{? > 0,M > 0} and the equihbrium (ii,Mi) = (0.5,0) is a saddle point. The 
third equilibrium point has U2 < and has no biological meaning. 

The above examples suggest that the equilibrium (ii,ui) is always a saddle 
and that (i2, U2) is locally asymptotically stable, if U2 > 0. If (12, 1*2) exists and 
is locally asymptotically stable, then its basin of attraction is larger for larger 
values of the transmission rate t (this is the intersection of regions (B) and (C) 
in figure [1]). 

3 Infection with multiple strains 

Assume that two Wolbachia strains A and B are present in the population and 
let lA, is and iab denote the number of individuals singly or doubly infected, 
in addition to u, the number of uninfected individuals. We use the same scaling 
as in the previous section, where time was scaled with respect to the birth rate 
h. There are seven incompatibility levels qx.Y where X,Y ^ {Q,A,B,AB} are 
the infection statuses in the incompatible cross. These are shown in the directed 
graph of incompatibilities in figure [3] 

We assume that the transmission of one strain in doubly infected individuals 
is independent of the transmission of the other strain ,23^. Moreover, the mor- 
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Figure 3: The directed graph of possible incompatibihty relations. An arrow 
from node X to node Y indicates that a Xd' x F5 cross is incompatible, with 
incompatibility level qY,x- 



talities due to infection with either strain are additive. With these assumptions, 
and setting p = iab + ia + is + u for the total population, our model is 



TATBlAB - VABPtAB, 



dtAB 

dt 

dlA \- s (^ Mb\ . 

-7- = Ta{1 - TB)tAB + TA 1 - qA,B qA,AB «A - VAPlA, 

dt \ P P J 

dlB ,^ \ • , /i iAB\ . 

—J- = (1 - TAjTBlAB + Ts 1 - qB,A qB,AB «B - VBPlB, 

dU ,^ N/-, ^ . , A / -, *S lAB ^ 



+ (1 


-rA}{ 


f 


iA 




iAB 




qB,AB 




P 


P 


IB 




iAB 


3 


- qo,AB 




P 


P 



J, \' ■^/\' -aj-i^a . V- i^,B qA.AB «A 

at \ P P 

%B 

P P J 

, , , I A T-B «Ab\ 

+ ( 1 - qo^A — - qo,B — - qo,AB ]u-pu, 
where 

d + DA d + De , d + DA + De 

Va = — , Vb = J and tjab = , = r?^ + 775 - 1 

d d d 

are measures of the fitness costs of the individual infection types. This model 
can now be reduced in its complexity in a variety of ways. For example, if the 
initial condition lies in the space {0} x R^, that is, there are no doubly infected 
individuals present initially, then the solution will also lie in that space at all 
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times. Moreover, setting appropriate incompatibility levels qx,Y to zero allows 
to study cases of mutual compatibility. 

As a first illustration, we want to consider the absence of doubly infected indi- 
viduals, iAB = 0, mutual compatibility of infected individuals, qA,B = qs.A = 0, 
equal transmission efficacy ta = tb ='■ t and equal infection costs rjA = Vb =■ V 
(again, we write £, = i]^^)- The symmetry is only broken by the levels of incom- 



patibility go, A 7^ qo.B- Hence, equation (10) simplifies to 

diA I . . 

diB , ^. , , 

-j^ ^ (T -r]p)iB, (11) 

du ( 14 iB\ 

-J- = {1-t){ia + ib) + [I- qa.A- qo.B — ] u-pu. 

dt \ P P J 

Observe that model ( [IT| ) has the property that planes orthogonal to the (m, is)- 
plane 

Ra^ { {iAi iB,u) e M>o : U - ais = 0} 



are invariant under the flow generated by (11). This is seen from the fact that 
for every a e [0, oo] 

^{iA - = (t - 'np){iA - aiB) = 

on Ra- This implies that the ratio ^^'^^^ remains constant along a trajectory. In 
other words, if transmission efficacies and death rates are equal for two strains 
(as are birth rates throughout in our model) then neither strain can replace 
the other among the infected individuals based on stronger cytoplasmic incom- 
patibility. This is in line with recent theoretical predictions of Turelli ,22j and 
Haygood and Turelli >9,, who suggest that strains are selected for relative fe- 
cundity rather than high levels of cytoplasmic incompatibility. It needs to be 
pointed out however, that even a small difference in transmission efficacies or 
death rates of the two strains helps the strain with the greater transmission rate 
or the lower mortality to establish itself in the population. 



The disease free equilibrium of ( 11 ) is easily found to be 



(«A,o, is.o, wo) = (0,0,1). (12) 

It is clear that the subspaces {0} x M x M and M x {0} x M are forward invariant 
under the flow generated by (111 and that the equilibrium solutions ([t]) exist in 



the respective subspaces, with g in ([t]) replaced by either go, A or qo,B- That is 
we have equilibrium solutions 

(iA,i(9o,A,T,0,0,T^ - iA,i), (iA,2('7o,A,T,^),0,T5 - ^,2), 
(0, iB,i{qo,B,T, 0, - «B,i), (is,2(0, qo,B,T, C), - »S,2)- 
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Besides that, it can be checked that there is a continuum of equilibrium solutions 



iB{u) 



— ')';Q,BTtl + 7j^go,Btt^+T^(l + (t7— 1)m) 
'^^(90,^1-90,3)" 

^riqo,Aru-ri^qo^Ay-^-T^{l+{ri-l) 
ri^iqo,A-qo,B)u 



(13) 



for every u G (0, r^), provided that these expressions are non-negative, 
solutions from ( 13 1 satisfy 



The 



iA{u) + iBiu) = U = — u, 



which corresponds to equation Q. 

Example 3.1. We consider the case r = 1 of complete transmission. One checks 



by direct calculation that system ( 11 ) then has another manifold of equilibrium 
solutions 

M + is = u = 0. 

The intersection of this line with each plane orthogonal to the (ij^, is)-plane 
is an equilibrium for the flow restricted to that plane. In addition, each Ra 
contains a saddle point. A numerical example is shown in figure |4] (left panel). 
If, in contrast, different costs are associated with the infection, the strain with 
the lower cost will dominate the population, see figure [4] (right panel). Similarly, 
if infection costs are equal, but one strain transmits more efhciently, then it is 
going to dominate the population. 



We return to model ( 10 ) and consider mutual incompatibility of singly in- 



fected individuals. The equations are 



iA + iB\ 

1 - (7o ]u-pu. 

P J 



1 - IB. A— «B 
P 



(14) 



It is clear that if either strain is not present initially then it will remain absent 
at all times. On the marginal spaces {ia — 0} and {ib = 0}, equations (14 1 
reduce to the single strain model ([l])-([2]) and have the corresponding equilib- 
rium solutions where only one strain is present, with the common disease free 
equilibrium (12 1. It follows from (fTIl) that 
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Figure 4: (Left) Dynamics of model (111 in the case that go, A = 0.95, go,B = 
0.5, T ~ 1 and ry = 1.1. The soUd blue line is a family of attractors for the 
flow restricted to planes Ra, the red line is a branch of saddle points in each 
of these subspaces. Several individual trajectories are shown, and the space 
Ri — {"lA = *b} is marked. (Right) If we choose instead r/A = 1-1 and rjB = 1.2 
(while keeping all other parameters the same) then the less costly strain A 
dominates the population. 



^{qB,AiA - qA.Bis) = qs.A - 1A.B^~ ] - VP ] ^A 

- qA,B (^T ^1 - qB,A^-^^ - VP^ IB 

= T (^qB.A (^1 - qAM^-^^ iA - qA,B (^1 - ^B.A^-^ I «B 

- VP{qB,A-iA ~ qA,B-iB) 

= (t - r]p) {qB,AiA - QA.BiB)- 

This implies that the plane orthogonal to the (i^i, iB)-pla'iie, 

-R* = {(«A,»B,m) e K>o '■ qB.AiA ^ QA.BiB} , 

is forward invariant, and hence so are the wedges on either side. Solving equa- 



tions ( 14 1 for the total population yields that if ia ^ 0, ib ^ 



= p p+ -qA.BiB, 

= p p+ -qB.AiA- 

rj ri 
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For this system to be consistent, it is necessary that if 7^ 0, is 7^ then 



qs.A^A = qA,BiB, 

hence any coexistence equilibrium of the two infected strains has to lie in the 
plane i?* . Indeed, there may be a saddle point equilibrium solution in i?^ that 
is locally asymptotically stable for flows that start in i?*. 

Example 3.2. Let qA,B = 0.99, qB,A and go, A = 9o,s = 1- Choose the trans- 
mission efhcacy t — 1 and the cost of the infection 77 = 1.1. We see in figure [5] 
that every solution starting outside the space i?* converges to an equilibrium in 
one of the marginal spaces. The plane i?* contains a locally stable equilibrium 
for trajectories starting in i?* which has M2 = and it contains a saddle point for 



trajectories starting within i?* (not shown, compare to example 2.1 and figure 
[21). 




Figure 5: Mutually incompatible strains as described in (14) do not show co- 
existence. The marginal equilibrium solutions are marked by dots, the disease 
free equilibrium (equation (12), green) attracts some solutions from the plane 
i?* and other trajectories beginning in the wedges on either side that are not 
shown here. 



Again, we need to recall that this dynamical behaviour is not generic, in 
the sense that the complement of the set {ta — tb, fjA = '7b} is dense in the 
parameter space. 



Finally, we want to explore the full model ( 10 ) when doubly infected indi 



viduals are present. To reduce the number of parameters somewhat, we assume 

qA,B = qA,AB — ?0,S, qB.A — qB.AB — qo,A, 

that is, the presence of one strain in the fertilizing male that is missing in the 
female has the same effect, regardless of the other infections that the female may 
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carry [101 P-66]. Moreover, the escape from cytoplasmic incompatibility for the 
offspring of an uninfected female and a doubly infected male is the product of 
the two individual escape probabilities [25 , 

1 - qo.AB = (1 - go,A)(i - qq^b)- 

Example 3.3. Under the above assumptions, let Qa.b = 0.9 = qB,A, the trans- 
mission efficacy ta — tb — 0.9 and the cost of the infection tja = tjb ^ l-l- 
Then there exists a coexistence equilibrium of doubly infected and both types 
of singly infected individuals, where however the proportion of doubly infected 
individuals is much larger (figure [6]). 



IB 



0.5 / 




0.5 



1 

0.5 



u 



Jab 



0.5 



Figure 6: With the parameters chosen as described in example 3.3 we observe 
stable coexistence of all three infection types. Solid, dashed and dotted lines, 
respectively, show parts of the solution starting at the same initial values. 



4 Introduction of age-structure 

In the previous sections we have seen that infection with Wolbachia gives rise 
to interesting dynamic behaviour already in unstructured populations. Clearly, 
individuals of different ages are subject to different fertility and mortality rates. 
We therefore expand our model ([l])-([2]) to include age-dependent fertility and 
mortality rates for infected and uninfected individuals. This leads to nonlinear 
partial differential equations with nonlocal boundary conditions that represent 
the birth process [5H81I24]. Although this results in more complex models, they 
are still amenable to analytical study. Here we focus on qualitative questions, 
when analytical progress is possible; in particular, how do the stability results 
for equilibrium solutions compare to the unstructured case. 

Let i{a, t) and u{a, t) denote the densities of infected and uninfected individ- 
uals of age a at time t, respectively, where a £ [0, m] (this is not to be confused 
with the notation in section |2j where they denoted scaled numbers of infected 
and uninfected individuals). Then the evolution of the population is governed 
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by 



itia,t)+iaia,t) = -r]i{a){I{t) + U(t))iia,t), 
ut{a, t) + Uaia, t) = ~T]2{a){I{t) + U(t))u{a, t). 



i{0,t) — T / I3i{a)i{a,t) da, 



(15) 
(16) 

(17) 



u(0,t) = (1 - t) / /3i{a)i{a,t)da 

m 



1-9 



I{t) + Uit)J J, 



I32ia)u{a,t)da, (18) 



where 



/(t) 



i{a,t)da, U{t) 



u{a, t) da, 



and r}i, ri2, Pi and ^2 denote the age-specific mortahty and fertihty rates for 
infected and uninfected individuals, respectively. System (15|-(18) is equipped 
with initial conditions 

i(a,Q) — ii^{a), u{a,0) — uo{a). 
The parameters r and q have the same meaning as in section [2j 

4.1 Existence of equilibrium solutions 



We find the time-independent solutions of equations (15 1 and (161 as 



i*(a) = i*(0) cxp | — (/* + ?7*) J rii{r) dr 
w* (a) = u* (0) exp <—(/* + t/* ) / rj2{r)dr 



+ 1-9- 



(19) 
(20) 

(21) 



u,{0) I /32 (a) exp <j + / 772 (r ) dr da. 

(22) 



where z*(0) and li^(O) satisfy 
i*(0) =Ti,(0) / /3i(a)exp<j -(/, +C/,) / m{r)dr}da, 



u,{0) = (l - T%iO) A(a)exp<^ -(/,+[/,) / r]i{r)dr}da 



Here 



(a) da, 



u*(a) da. 
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First we note that the trivial steady state (0, 0) always exists. Next we note 
that if = then equation (22) reduces to 



1=1 /32(a) exp I -t/* J rj2{r)dr^ da. 



(23) 



It then immediately follows from the monotonicity and continuity of the right 



hand side of (23) (as a function of [/*) and the Intermediate Value Theorem, 
that a unique disease free equilibrium, given by 



M*(a) = 



f/* exp{-J7* /o" ??2(?') dr} 
/™ exp {-[/* 7/2 (r) dr} da ' 



exists if and only 



^2(a)da > 1 



(24) 



holds. If we look for strictly positive equilibrium solutions (i*(a), u*(a)) we find 
that and C/* have to satisfy 



(25) 



l=Ty^ /3i(a) exp I -(/, + [/,) 7?i(r)dr|da, 

U* (1 - (1 - Qt:^) r /32(a) exp {-(/* + U.) m{r) dr) da) 
rexp{-(/,+C/,)/;ry2(r)dr}da 
/* (r-i - 1) 



/o" exp { - (/* + ) ?7i (r) dr } da ' 



(26) 



Conversely, if and satisfy equations (p5|)-(26) then equations (19)-(20l 



determine uniquely a positive equilibrium solution. We also see from equation 



(25) that 



(3i{a)da > — , 



(27) 



is a necessary condition for the existence of a positive equilibrium. In fact, if 



equation (271 holds true, then we can solve equation (251 to obtain a unique 



positive value 



(28) 



A straightforward calculation then leads from equation (26) to the following 
quadratic equation for 



I^Cs +/*(!- C2 - C1C3 + C4) + C1C2 - Ci = 0, 



(29) 



where we have defined 

C2 = / /32(a) exp ■{ -ci / r/2(r) dr J> da. 











C3 



qc2 

Cl ' 



and 



C4 = 



^ ~ ^) /o" {-^1 lo '72(?')dr} da 
/o" <5xp {-Cl /° r]i (r)dr} da 
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Similarly to the unstructured case, see equation ([5]), we arrive at a quadratic 
equation (unless g = 0) for the infected population size /*. Of course the cal- 
culations now are much more involved since we have age-dependent fertility 
and mortality rates. However, for fixed model ingredients the equilibrium so- 
lutions can be determined explicitly, via equations (19)-(20). In contrast to 



the unstructured case, we have necessary conditions on the birth rates for the 
existence of non-trivial equilibria. 

We summarize our findings in the following theorem. 



Theorem 4.1. The equilibrium solutions to equation system ( 15 I -( 18 1 are given 



by functions (19) -(20) with initial values (21) -(22), provided that the total pop 



Illations of infected and uninfected individuals and given by equations ( 29 I 



and (28) are non-negative. 



Wc note the formal similarity of equations ( 28 1 and ( 29 ) to the conditions 
([4| and ^ for the unstructured model in section [2[ 

4.2 (In)stability 

In the previous section we established necessary and sufficient conditions for the 



existence of a non-trivial steady-state of system (15)-(18|. In this section we 



study stability properties of the steady states. To this end, first we formally 
linearise system (15)-(18) around a steady-state solution (i*(a),u*(a)). We 
introduce the perturbations p{a, t) = i{a, t) — i^,{a) and s(a, t) = u{a, t) — Ut(a) 
and we use Taylor series expansions of the fertility and mortality functions. 
Then we drop the nonlinear terms to arrive at the linearised system 



Pt{a,t) +pa{a,t) 
St{a,t) +Sa{a,t) 



-m{aKp{a,t){h + U, 
-r]2{a){s{a,t){I.^. + 



PiO,t) 
s{0,t) 



l3i{a)p{a, t) da. 



Jo 



Pi {a)p{a, t)da + [1 — q 



^ p(t] 



-iM{Pit) + sm, (30) 

u,{a)iP{t) + Sm, (31) 

(32) 

/?2(a)s(a, t) da 

(33) 



Sit) 



U* I Jo 

/32(a)ii*(a) da 



where 



P{t)= / p(a,i)da, S{t) = I s{a,t) da. 
Jo Jo 



For more detailed calculations we refer the reader to [SHZ], where similar type 
of age and size-structured models where treated. It can be shown that the 
linearised system is governed by a strongly continuous semigroup of linear op- 
erators, which is eventually compact (see e.g. [HIIZ])- However, this governing 
semigroup cannot be shown to be positive, since mortality of both infected and 
uninfected individuals is an increasing function of the total population size. 
Eventual compactness of the governing linear semigroup implies that to study 
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stability of steady-states it is sufficient to study the point spectrum of the lin- 
earised operator see, e.g. The standard way how this can be carried out is to 
solve the eigenvalue equation and deduce a characteristic equation (if possible) 
whose roots are the eigenvalues of the linearised operator. We note that the 
lack of positivity implies that we cannot expect to establish even local stabil- 
ity results unless the characteristic equation can be cast in a simple form. We 
substitute 



p{a,t) 
s{a, t) 



exp {At} 



v{a) 
w{a) 



into the linearised equations p0|-(33). This yields 



v'{a) 
w'{a) 

«(0) 
w(0) 



-v{a) (A + mia){I* + U,)) - m{a)iM{V + W), 
-w{a) (A + m{a){h + U,))- m{a)u,{a){V + W), 



l3i{a)v{a) da. 



(1-r) I3i{a)v{a)da+ (^1 - g-^ 
1~F7 — TT \9 / P2[a)u^[a)da, 



(34) 
(35) 

(36) 



P2{a)w{a) da 



where 



V = 



(37) 







v{a) da, 



I w{a)da. 
/o 



Hence A € C is an eigenvalue if and only if the nonlocal system (|34|)-(p57|) admits 
solution. The solutic 

via)^fl{a) (viO) 



a non-trivial solution. The solution of the differential equations (34) and (35) is 

7ji{x)i.,{x){V + W) 



w{a) = flia) ( w{0) 
where we have introduced 



fli^) 
" r^(x)w,(a;)(F + M^) 





dx 
dx 



(38) 
(39) 



/1(a) = exp - / A + n,{y){h + C/*) dy 



1,2. 



Next we substitute the solutions (38) and (39) into the boundary conditions 



( 36 ) and ( 37 1 and integrate the solution ( 38 ) and ( 39 ) from zero to m to arrive 



at a four-dimensional homogeneous system for the unknowns v(^),w(fS),'V and 
W . This homogeneous system admits a non-trivial solution if and only if the 
determinant of the coefficient matrix equals zero. We can formulate the following 
theorem. 
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Theorem 4.2. X is an eigenvalue of the linearised operator if and only if it 
satisfies the equation 



I TflsCA) - 1 



K{\) = det 







-rag (A) 



(l-r)a5(A) ariX)~^^-l a^iX) 



aiiX) 
\ 





asiX) 



-TaeiX) \ 
aio(A) 

-a2(A)-l -a2(A) 
-a4(A) -a4(A) - 1 / 



where 
ai{X) 

asiX) 

a^iX) 

ajiX) 



(40) 



/a (a) da, 
/a(«) da, 



02 (A) 
a4(A) 



/a(«) 
/a(«) 



" ?7i(a;)i4a^) 
" 772(2;)"* (a;) 



da; da, 
dx da, 



/3i(a)/A(a)da, aelA) 
I32{a)fl{a)da, aQ{X) 



Pi{a)flia) 
' P2{a)fl{a) 



° T]i{x)i.^{x) 
•2 

A 



Ax da, 
dx da, 



(A) = (r - l)ae(A) + - 1 j a.(A) - ^^-^ ^.(a)..(a) 



) da. 



{h + u^Y Jo 



aio(A) = (r - l)a6(A) + [^jff^^ - 1 j a9(A) + 
It follows from the growth behaviour of the functions f{ that 



/32(a)M*(a) da. 



lim K{X) = det 

A— >+oo 



/-I \ 

-1 Ci C2 

0-10 

yo 0-1/ 



(41) 



the limit being taken in M, and Ci,C2 are constants. Hence we can formulate 
the following general instability criterion, which follows immediately from the 
Intermediate Value Theorem. 



Theorem 4.3. The stationary solution (i^(a), u^(a)) of equations (15)-(18) is 
unstable if K{0) < 0. 



As we can see the characteristic equation (40) is rather complicated in ge- 
neral, hence we only consider some interesting special cases when analytical 
progress is possible. 
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4.2.1 The trivial steady state 

We consider the stability of the steady state i^, = 0, li* = 0. Note that in this 
case the characteristic equation (|40|) reduces to 



K{X) = dct 



/ rar^iX) - 1 \ 

(l-T)a5(A) a7(A)-l 

ai(A) -10 

V asiX) -1 J 



= 0, (42) 



which leads to the equation 

{Ta,{X) - l)(a7(A) - 1) = 0. (43) 

Therefore, A G C is an eigenvalue if and only if A satisfies either of the two 
equations 

l = T /3i(a)e^^''da, 1= P2{a)e-^'' da. (44) 
^0 Jo 

We can formulate the following theorem. 

Theorem 4.4. The trivial steady state is locally asymptotically stable if 



/3i(a)da<l and / /32(a) da <1. (45) 
Jo 

On the other hand, if either 

I pTn 

/3i(a)da>l or / /32(a)da>l. (46) 

"'0 

holds, then the trivial steady state is unstable. 
4.2.2 The disease free steady state 

Consider the disease free steady state, i.e. = 0, which exists by condition 
(24) if and only if J™ /32(a) da > 1. In this case the characteristic equation (40 1 
can be written as 

K{X) - (ra^iX) - 1) 

a7(A)-l _ iiT^^^lliiM^ __a9(A) 



(47) 



X det I -1 

a3(A) -a4(A) -a4(A) - 1 

= (rasiX) - 1) [(-a4(A) - l)(a7(A) - 1) + a3(A)a9(A)] 
= 0. 

This again splits into two equations. The first one is easy to analyse, since it 
can be written as 



l = r^ /3i(a)exp|-;7, t]i{x) dx^ e-^" da. 



(48) 
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Theorem 4.5. // 



/3i(a) exp < — [/* / 771(2;) da; > da > 1 



(49) 



where U* satisfies (23 I, then the disease free steady state is unstable. 



Remark 4.6. Provided that equation (47) has a dominant real solution X, it is 
shown that condition (49) in Theorem 4-5 is indeed necessary and sufficient for 



the instability of the disease free steady state. However, as we have noted before 
the governing linear semigroup cannot shown to be positive, hence we cannot 
establish the existence of a dominant real root of the characteristic function 



(40) 



4.2.3 Complete transmission of disease 

In case of complete transmission of disease, i.e. when r 
can be written in the following form 



1, equations ( 15 )-( 18 1 



where 



it{a,t) + ia(a,t) 
ut{a,t) + Ua{a,t) 

i{0,t) 
u{0,t) 

fi,{a,I{t),Uit)) 
/3|(a,/(t),C/(i)) 



-fiiia,I{t),U{t))i{a,t), 
-fi2ia,I{t),U{t))u{a,t), 

l3i{a)i{a, t) da, 
Pl{a,I{t),U{t))u{a,t) da. 



T1,{a){I{t) + U{t)), i = 

m 
'm + u{t) 



1,2, 



/92(a) [l-q 



(50) 
(51) 

(52) 
(53) 

(54) 
(55) 



Hence model (|50[)-( 53 ) is a special case of the n-species age-structured system 
considered in [ojj where the coupling occurs due to competition for resources 
and due to the inhibition of the proliferation of the uninfected population. 

In [S] we deduced a very general instability condition, which we recall for 
the case n = 2 for the readers convenience (see Theorem 2.3 in [5])- 

Theorem 4.7. A strictly positive .stationary solution (i^,u^,) G (I 



'>0j 



of (50) 



( 53 ) is unstable if the partial derivatives of the net reproduction rates of the 

(56) 



infected, respectively uninfected populations satisfy 



Taking into account (54)-(55) we have 



R^{I,U)=J^ /3i(a)exp|-(/ + f7)^ 77i(x)da;| da, (57) 
R\l,U)^j^ /?2(a) (^l-(7-^^exp|-(/ + C/)^ 772(x)da;| da. (58) 
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From equation ( 57 1 wc obtain 



r]i{x) dx ) exp <—(/*+[/*) / 771(2;) dx > da < 



unless /3i = or 771 = 0. Also from equation ( 58 ) we obtain 



/32(a) [l-q 



j T]2 (x) dx^ exp I - (/* 



/32(a) exp <j -(/, + [/,) / 772(3;) da; )■ q 



772(2;) dx > da 



da, 



/32(a) \l~q 



j ri2 (x) dx^ exp | - (/* 



+ 1 /32(a) exp <{-(/,+[/*) / 772 (a;)da; U 



772(2;) dx > da 



(/*+[/*)■ 



da. 



Hence Rfj{I^,U^) > R^j{h,U^,) for every strictly positive steady state, unless 
/32 = or g = (in which case is the strictly dominant eigenvalue of the 
linearised operator). We summarize our findings in the following theorem. 

Theorem 4.8. Assume that t — 1, q ^ Q, and /3i,/32,?7i are not identically 
zero. Then any strictly positive steady state of equations (15) -(18) is unstable. 

In other words, there is no coexistence of infected and uninfected popula- 
tions, This corresponds to the instability of the equilibrium solution (ii,iti) in 
the left panel of figure [2] for the unstructured case. 



5 Discussion 

In the present work we introduced and studied differential equation models for 
the dynamics of populations infected with Wolbachia. First we built ordinary 
differential equation models, in which we have implemented fitness costs of an 
infection as increased mortalities while keeping the birth rates equal for all 
infection statuses. It is equally appropriate to reduce birth rates for infected 
individuals and (for the sake of simplicity) then to assign the same mortality 
to all individuals. This leads to the following model for the case of a single 
Wolbachia in an asexual population 

di 

^ = (/^T- {i + u))i, 

^ = ^(1 - r)i + ( 1 - q—^ ii + u)] u, 

at \ I + u J 
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where /i G [0, 1] is the reduced fecundity of infected individuals. This results in 
similar formulas for equilibrium solutions as ([t]) and vector fields as in figure [2j 
An experimentally testable prediction of our model ([T])-([2| is that there are no 
persistent Wolbachia strains with a transmission efficacy less than | (see the 
region for existence of the observed stable equilibrium (^21*2) in figure [T]). 

Our model for multiple infections is novel insofar it allows the theoretical 
biologist to adapt it to a case of special interest (with or without doubly in- 
fected individuals, with or without mutual incompatibility). This should help 
to gain a more unified perspective than what was possible from models cre- 
ated for each purpose individually. In the case of mutual compatibility we saw 
that strains with higher transmission efficacy or lower mortality due to infec- 
tion establish themselves over competitors. This is in good agreement with 
other predictions from discrete population genetics models |9l [22]. Although 
the model for infections with multiple mutually incompatible strains in section 
[3] is too complicated for all of its equilibrium solutions to be written down ex- 
plicitly, it can be analyzed to a certain degree by identifying invariant subspaces. 
By numerical simulations we provided evidence for the absence of coexistence 
of singly infected individuals, apart for exceptional choices of parameters and 
initial values. This situation changes, if doubly infected individuals are present 
that can lose one of their strains when giving birth. 

We have expanded the simple ordinary differential equation model from sec- 
tion[2]by introducing age-structure. Of the extensive literature about structured 
populations, let us mention the monographs [2 [13 [23], the classical paper [S] 
and the recent collection [13] . Clearly, age-structured models allow a much finer 
level of detail to be incorporated, but also pose greater analytical challenges. 
Nevertheless, we have shown that existence of equilibrium solutions and their 
stability properties can be investigated in a straightforward fashion. We saw 
that unlike in the unstructured case, existence of a disease free equilibrium is 
now subject to a condition on the integral of the birth rate. We also obtained 
analytical results in some special cases which allow at least a partial character- 
ization of the dynamic behaviour of the system once the model parameters are 
fixed. 

Wolbachia together with Cardinium are the two bacterial infections of arthro- 
pods that cause cytoplasmic incompatibility. In this work we have focused on 
the case of diplodiploid species where cytoplasmic incompatibility results, with 
a certain probability, in embryonic death. In future work, we will include sep- 
arate sexes into our models, see the book by lannelli et al. [11] for a compre- 
hensive introduction to gender-structured populations. Then it will be possible 
to study gender-specific effects of the Wolbachia infection. These become even 
more important in haplodiploid organisms such as bees, ants and wasps, where 
cytoplasmic incompatibility is vastly more complex [19|,l23j (male-development, 
thelytokous parthenogenesis etc.). 
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